/*
 * Voronoi_channels.cpp
 *
 *  Created on: Jun 8, 2012
 *      Author: marchi
 */

#include "Voronoi_channels.h"

namespace Cavities {
double isLargerR::radius=0.0;
void Voronoi_channels::getDelaunay(isTypes & isWhat){
	int ia=0;
	if(!Mycon) return;
	c_loop_order_periodic vl(*Mycon,*porder);
	voronoicell_neighbor c;
	vector<double> vertx;
	cv.clear();
	cvnei.clear();
	cvnd.clear();
	AtRadius.clear();
	set<Vertex> cvE;
	if(Atoms.empty())
		Atoms=vector<Dvect>(nr,0.0);

	if(vl.start())
		do {
			vector<int> nodes;
			if(Mycon->compute_cell(c,vl)) {
				int p=c.p;

				c.vertices(vertx);
				double x,y,z,r;
				int iD;
				vl.pos(iD,x,y,z,r);
				Atoms[cindex[ia]](x,y,z);
				AtRadius.push_back(r);
				vector<Vertex> v0;
				vector<int> iv0;
				for(int m=0;m<p;m++){
					int xx=m*DIM+XX;
					int yy=xx+1;
					int zz=yy+1;
					double myr=sqrt(vertx[xx]*vertx[xx]+vertx[yy]*vertx[yy]+vertx[zz]*vertx[zz]);
					nodes.push_back(cindex[ia]);
					if(c.nu[m]>0) {
						int j=0;
						while(j<c.nu[m]-1) nodes.push_back(cindex[c.ne[m][j++]]);
						nodes.push_back(cindex[c.ne[m][j]]);
					};

					try{if(nodes.size() != 4) throw string(" Warning Number of nodes is not 4 but ");}
					catch(const string & s){cout << s << nodes.size() << endl;}

					Vertex tm(nodes);
					Vertex tm2(tm,Top.gbG());
					if(myr-r < isLargerR::radius) {
						cvE.insert(tm);
						map<const Vertex,CavityData >::iterator it=cv.find(tm);
						if(it != cv.end()) cv.erase(it);
					} else if(cvE.find(tm) == cvE.end() && isWhat(tm2)){
						cv[tm](myr-r,Dvect(vertx[xx]+x,vertx[yy]+y,vertx[zz]+z),tm2);
						iv0.push_back(m);
					}
					v0.push_back(tm);
					nodes.clear();
				}
				for(unsigned int ma=0;ma<iv0.size();ma++){
					int m=iv0[ma];
					for(int j=0;j<c.nu[m];j++) {
						Vertex tm=v0[c.ed[m][j]];
						if(cvE.find(tm) == cvE.end())
							cvnei[v0[m]].insert(tm);
					}
				}

			}
			ia++;
		} while(vl.inc());
}
double Voronoi_channels::intersect(Dvect & la,Dvect & lb,vector<int> & nei){
	Dvect p0=Atoms[nei[0]];
	Dvect p1=Atoms[nei[1]];
	Dvect p2=Atoms[nei[2]];
	double r0=AtRadius[nei[0]];
	double r1=AtRadius[nei[1]];
	double r2=AtRadius[nei[2]];
	Dvect xab,x10,x20,xa0;
	xab=la.Minus(lb,co,oc);
	x10=p1.Minus(p0,co,oc);
	x20=p2.Minus(p0,co,oc);
	xa0=la.Minus(p0,co,oc);
	Matrix mx;
	for(unsigned int n=0;n<DIM;n++){
		mx[n][XX]=xab[n];
		mx[n][YY]=x10[n];
		mx[n][ZZ]=x20[n];
	}
	Matrix mix=mx.Inversion();
	Dvect t=mix*xa0;
	Dvect pos;
	for(unsigned int n=0;n<DIM;n++){
		pos[n]=la[n]-xab[n]*t[n];
	}
	set<double> rda;
	rda.insert(pos.Dist(p0)-r0);
	rda.insert(pos.Dist(p1)-r1);
	rda.insert(pos.Dist(p2)-r2);
	return *rda.begin();
}
bool Voronoi_channels::Percolate(mapVCD::iterator & itc, mapVCD & v, mapVCD & cv){
	mapVsV::iterator itver=cvnei.find(itc->first);
	v.insert(*itc);
	Dvect la=itc->second.gPoint();
	vector<int> nei_a=itc->first.gnoVertex();
	cv.erase(itc);
	if(itver == cvnei.end())return false;
	set<Vertex>::iterator itsver=itver->second.begin();
	for(;itsver != itver->second.end();++itsver){
		mapVCD::iterator itcnew=cv.find(*itsver);
		if(itcnew != cv.end()) {
			Dvect lb=itcnew->second.gPoint();
			vector<int> nei_b=itcnew->first.gnoVertex();
			vector<int> mynei=pickPlane()(nei_a,nei_b);
			double r=intersect(la,lb,mynei);

			if(r >= isLargerR::radius){
				return Percolate(itcnew,v,cv);
			}
			else
				return true;
		}
	}
	return true;
}

void Voronoi_channels::Channels(){
	map<Vertex,CavityData > newcv;
	map<Vertex,CavityData >::iterator it=cv.begin();
	for(;it != cv.end();)
		if(isLargerR()(*it))
			newcv.insert(*it++);
		else ++it;


	it=newcv.begin();
	mapVCD vv;
	mChannels.clear();
	while(it != newcv.end()){
		if(Percolate(it,vv,newcv)){
			mapVCD::iterator itb=vv.begin();
			bool bW0=false;
			bool bW1=false;
			for(;itb != vv.end(); ++itb){
				vector<int> tm=itb->second.gVtx().gnoVertex();
				if(count_if(tm.begin(),tm.end(),bind1st(equal_to<int>(),0))) bW1=true;
				else bW0=true;
			}
			if( bW0 && bW1){
				mChannels.push_back(vv);
			}
			vv.clear();
		}
		it=newcv.begin();
	}
}

void Voronoi_channels::WriteIt(std::ofstream & fout){
	vector<mapVCD>::iterator it=mChannels.begin();
	for(unsigned int n=0;n<mChannels.size();n++)

	vector<vector<int> > tmp=Top.gvG();
	for(int n=0;n< Top.gbGMax();n++){
		double sum=hstg[n].TotHist();
		hstg[n].Avg();
		string ss=" ";
		ss+=Top.AtomResidue(tmp[n][0])+ " : ";
		if(n)
			for(unsigned int m=0;m<tmp[n].size();m++){
				ss+=label[tmp[n][m]]+ " ";
			}
		fout <<"& " << n << "  sum = " << sum << ss << endl;
		fout << hstg[n];
	}
}

} /* namespace Cavities */
